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SUMMARY 


A  finite  difference  code  for  predicting  the  high-speed  flow  over  the 
advancing  helicopter  rotor  is  presented.  The  code  solves  the  low-frequency, 
transonic  small  disturbance  equation  and  is  suitable  for  modeling  the  effects 
of  advancing  blade  unsteadiness  on  blades  of  nearly  arbitrary  planform.  The 
method  employs  a  quasi-conservative  mixed  differencing  scheme  and  solves  the 
resulting  difference  equations  by  an  alternating  direction  scheme.  Computed 
results  show  good  agreement  with  experimental  blade  pressure  data  and  illus¬ 
trate  some  of  the  effects  of  varying  the  rotor  planform.  The  flow  unsteadiness 
is  shown  to  be  an  indispensible  part  of  a  transonic  solution.  It  is  also  shown 
that,  close  to  the  tip  at  high  advance  ratio,  cross-flow  effects  can  signifi¬ 
cantly  affect  the  solution. 


INTRODUCTION 


Air  flow  past  a  helicopter  rotor  blade  exhibits  many  very  complex  features 
such  as  three-dimensional  unsteady  effects,  shock-wave  motions,  vortex  inter¬ 
actions,  and  stall.  A  complete  numerical  simulation  cannot  even  be  attempted 
yet,  but  it  is  possible  with  the  present-day  computers  and  numerical  methods 
to  model  some  of  these  features  and  acquire  a  better  understanding  of  some  of 
the  mechanisms  involved. 

The  model  used  in  this  study  is  a  perfect  fluid  model  that  is  further 
simplified  by  the  small-disturbance  approximation.  Weak,  almost  normal  shock 
waves  are  accounted  for  by  retaining  the  leading  nonlinear  term  in  the  stream- 
wise  direction.  This  model  is  useful  for  simulating  the  subsonic  and  transonic 
flow  past  the  advancing  blade.  Under  these  conditions  the  incidence  is  usually 
small,  and  the  results  presented  correspond  to  nonlifting  blades.  A  proper 
wake  representation  is  required  to  extend  this  simulation  to  lifting  configur¬ 
ations.  Prediction  of  the  complicated  rotor  vortex  structure  is  not  within 
the  scope  of  the  present  work. 
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EQUATION  AND  BOUNDARY  CONDITIONS 


The  mathematical  model  used  in  this  report  is  the  three-dimensional 
unsteady  (low-frequency)  small-disturbance  transonic  equation  as  derived  by 
M.  P.  Isom  (ref.  1,  p.  20).  This  equation  is  derived  in  a  blade-attached 
Cartesian  coordinate  system  under  the  usual  assumptions: 

1  -  M2(l  +  p)2  -  0(62/3) 

e  *  0(6) 

tip  Mach  number  due  to  the  blade  rotation 
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inverse  of  the  aspect  ratio 
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In  condensed  notation  the  equation  can  be  written: 
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t,  x,  y,  and  z  are  the  dimensionless  dependent  variables  normalized  by  SI, 
c,  R,  and  6-1'3  c,  respectively,  and  y  is  the  ratio  of  the  specific  heats. 
At  each  time  step  NS  in  THREED,  the  coefficients  are  computed  and  stored  in 
one-dimensional  arrays 


A(J),  B(J),  BP(J),  C(J),  D(J),  E(J),  J  -  I,  .  .  .  JM 


for  all  values  of  the  spanwise  index  J.  Allowance  is  made  in  the  code  for  a 
term  A'  3<J)/3t  for  which  the  value  of  the  coefficient  is  stored  in  AP(J)  and 
has  been  set  to  zero  for  all  present  uses. 

Initial  and  boundary  conditions  are  required.  To  integrate  this  equation 
the  initial  condition  used  is  usually  the  quasi-steady  solution  (i.e. ,  <|>xt  ■  0 
in  eq.  (1)). 

On  the  mean  surface  of  the  blade  the  flow  tangency  condition  is  expressed 
(cf.  ref.  1)  as: 


■  (y  +  y  cos  t)f* (x)  at  z  ■  0 

At  the  innermost  grid  location,  ym^n»  two  boundary  conditions  can  be 
used: 

1.  A  symmetry  condition  (equivalent  to  a  flat  tunnel  wall  in  wing 
calculations) 


|^-0  (specified  by  setting  JSYM  -  1  in  THREED) 
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2.  A  strip-theory  condition  (used  for  rotor  blades  or  semi- inf inite  wings) 


3y2 
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In  the  far  field  a  Dirichlet  (4>  =  0)  or  a  Newman  (3<|>/3n  =  0)  condition  has  been 
used.  The  upstream  boundary  is  usually  taken  as  the  uniform  undisturbed  flow 
(♦  *  0). 


COORDINATE  TRANSFORMATION  AND  THE  CORRESPONDING  MESH  SYSTEM 


In  order  to  treat  a  large  class  of  planform  shapes,  a  coordinate  trans¬ 
formation  is  made  prior  to  the  discretization  of  the  equation.  This  transfor¬ 
mation  incorporates  some  one-dimensional  stretching  capabilities  concentrating 
the  mesh  in  regions  of  large  gradients;  in  particular,  near  the  surface  of  the 
blade,  near  the  leading  edge,  and  near  the  tip.  The  coordinate  transformation 
is  of  the  form,  (x,y,z)  -*  (5,n,S), 

where 


5  *  5(x,y) 


n  *  n(y) 
c  »  e(z) 


Equation  (1)  now  becomes: 
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(2) 

The  coefficients  in  equation  (2)  are  the  partial  derivatives  of  the  trans¬ 
formation.  This  form  of  the  equation  is  called  semiconservative;  the  metric 
coefficients  are  brought  outside  the  3/35,  3/3n»  3/3C  symbols.  It  can  be 
shown  that,  if  the  transformation  is  sufficiently  regular,  the  jump  conditions 
are  preserved  across  a  discontinuity. 

Computation  is  made  of  four  first  partial  derivatives  and  five  second 
partial  derivatives.  They  are  35/3x,  3£/3y,  3n/3y*  and  H/3z  (called  in 
THREED  XIX,  XIY,  YIY,  and  ZIZ,  respectively)  and  3*5/3x2,  325/3x3y,  325/3y2, 
32n/3y2,  and  32c/3z2  (called  in  THREED  XIX2,  XIXY,  XIY2,  YIY2,  and  ZIZ2, 


A 


respectively) .  These  quantities  are  computed  at  each  interior  mesh  point  by 
using  finite  difference  approximations  of  the  coefficients  of  the  inverse 
transformation  and  the  following  identities: 
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The  derivatives  3x/35»  .  .  .  are  evaluated  by  finite  differences  at  point 
i,j  as 


/3x\  Xi+i,j  ~  Xi-i,1 

2a? 


These  expressions  include  second-order  terms  [0(A5  +  An)2  +  0(A£)2].  The  mesh 
is  constructed  in  three  steps.  In  the  first  step  the  locations  of  the  span- 
wise  stations  are  defined.  The  following  stations  are  specified: 


ymin 


innermost  station  on  the  blade,  typically  y  -  0.5  (y  is 
referred  to  as  YN  in  THREED)  min  min 


y  a  special  station  on  the  blade  (e.g. ,  a  kink  in  the  planform) , 

c  typically  yc  =  0.9  (yc  -*■  YC  in  THREED) 

y^  the  tip  of  the  blade  yd  ■  1  (y^  -*■  YD  in  THREED) 
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y  outermost  radial  station,  typically  y  -  1.5  (  y  YX  in 

max  THREED)  max  x 


In  addition  to  these  real  numbers  the  corresponding  integers  (JC  <  JD)  must  be 
defined.  This  determines  how  many  stations  for  computation  are  located 
between  1  (ymin)  and  JC  (yc) ,  JC  (yc)  and  JD  (yd),  and  JD  (yd)  and  JM  (ymav) . 
The  following  analytical  expressions  are  used  to  define  the  mesh  stations: 


J  >  JD 


J  <  JC 

JC  <  J  <  JD 


Y (J)  =  YX  +  cos 


Y (J)  =  YN  +  cos 

n  -  n 


(5^) 


j(YD  -  YX) 


(YC  -  YN) 


Y (J)  =  YC  + 


n_  -  n 

D  c 


(YD  -  YC) 


where  the  variable  n  is  defined  between  0  and  1  by 

J  -  1 
n  “  JM  -  1 


The  planform  of  the  blade  then  yields  the  locations  xa  and  Xf  of  the  leading 
and  trailing  edges  as  functions  of  J.  For  this  purpose,  a  piecewise  analyti¬ 
cal  representation  is  made  of  the  planform.  The  chordwise  coordinate  trans¬ 
formation  has  no  radial  dependence  for  all  points  beyond  the  tip.  In  THREED 
xa  and  Xf  are  called  XA(J)  and  XF(J). 

The  second  step  in  mesh  construction,  in  the  chordwise  direction,  is 
defining 

Xn,in  upstream  boundary,  typically  x,^  -  -8  (x,,,-^  -*■  XN) 

^max  downstream  boundary,  typically  xmax  -  6  (x,^  ->  XX) 


and  the  indices  IA  <  IF  which  determine  how  many  stations  are  located 
between  1  (x,„fn)  and  IA  (xa) ,  IA  (xa)  and  IF  (xf ) ,  and  IF  (xf)  and  IM  (x,,,^) . 
Similar  analytical  expressions  are  used  to  define  the  mesh  stations  in  x: 

I  >  IF  X(I,J)  -  XX  +  cosfe  y-  €F)[XF(J)  -  XX] 


I  <  IA  X(I,J)  =  XN  +  cos/y  _j""^A  )[XA(J)  -  XN] 


IA  <  I  <  IF  X(I,J) 


XA(J)  + 


[XF(J)  -  XA(J) ] 
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where  the  variable  £  is  defined  between  -1  and  1  by 


C  =  -i 


,  2(1  -  1) 
IM  -  1 


In  the  third  step  in  the  vertical  direction  the  following  are  defined: 


zmin  lower  boundary,  typically  zmin  -  -3  (zmin  ->  ZN) 


zmax  uPPer  boundary,  typically  zmax  -  3  (zmax  ZX) 

and  the  indices  KU  =  KO  +  1  which  determine  how  many  stations  are  located 
between  1  (zmin)  and  KO  (nearest  to  the  lower  surface  of  the  blade),  and  KU 
(nearest  to  the  upper  surface  of  the  blade)  and  KM  (z^x)-  The  mesh  stations 
in  z  are  defined  by  using  the  analytical  expressions: 

K  >  KO  Z (K)  =  ZX  -  cos^-  c)zX 
K  <  KO  Z(K)  =  ZN  -  cos0  ^ZN 


where  z,  is  defined  between  -1  and  1  by 


.  2  (k  -  1) 
KM  -  1 


The  mesh  dimensions  in  the  code  are  set  up  to  allow  for  maximums  of  IM  -  64, 
JM  -  32,  and  KM  *  32. 


FINITE  DIFFERENCE  SCHEME 


In  equation  (1),  the  nonlinear  term  (3/3x)  [B(3<j>/3x)  +  B'  (3iJ>/3x)2] ,  which 
is  often  written  nonconservatively  as  V<J>XX,  is  responsible  for  the  mixed 
character  of  the  flow.  It  is  well  established  that  a  mixed  scheme  must  be 
used  for  the  nonlinear  flux  discretization  (refs.  2  and  3),  given  as  follows 
for  a  uniform  mesh  spacing: 


Let 


V.  =  B  +  2B' 
i 


i+i 


2  Ax 


i-i 


In  the  following  four  cases  to  be  considered  the  nonlinear  term  is  discretized 
e.g.  , 

Case  1  >  0  >  0  (subsonic  point) 


(the  indices 


j 


Discretization:  - 

and  k,  which  are  invariant. 


2*i  +  *i-i 
Ax2 

are  not  indicated) . 
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Case  2  <  0  <  0  (Supersonic  point) 


Discretization:  V . 


Pi  -  +  *i-2 


i-1  Ax2 

Case  3  <  0  -  0  (sonic  point) 


Discretization:  V. 


Pi  -  24>1_i  +  ^_2 

Ax2 


Case  4  >  0  V^_1  <  0  (shock  point) 


♦i+l  "  2<t‘i  +  $±-1  *±  “  2<t>i-l  +  l))i-2 

Discretization:  - r - —  +  V.  1  12 


Ax 


i-i 


Ax 


In  contrast  to  most  small  disturbance  codes  (typified  by  refs.  4  and  5), 
the  discretization  of  the  sonic  point  (case  3)  eliminates  some  spurious  oscil¬ 
lations  that  appear  when  the  sonic  line  is  located  close  to  the  leading  edge 
of  a  blunt  airfoil  in  a  region  where  the  flow  experiences  a  rapid  acceleration. 
It  can  be  shown  that  the  discretization  that  is  proposed  here  is  consistent 
with  the  equation,  but  it  is  not  strictly  conservative.  However,  the  error  of 
conservation  is  small,  and  not  larger  than  0(Ax).  The  shock-point  discre¬ 
tization,  however,  ensures  conservation  of  mass  at  the  shock  point. 

The  next  term  in  equation  (1)  is  the  cross-derivative  term.  This  term  is 
small  inboard  where  the  flow  is  subsonic  and  two-dimensional.  However,  for 
large  advance  ratios  (y  -  0.5)  and  for  values  of  azimuth  and  radius  where  the 
transonic  flow  has  a  large  radial  component,  its  effects  cannot  be  neglected. 

In  fact,  in  these  cases  the  cross-derivative  term,  which  is  usually  treated 
explicitly  (i.e.,  always  at  the  previous  time  level  (ref.  5)),  has  a  destabi¬ 
lizing  effect  and  can  strongly  reduce  the  time  step  required  for  maintaining 
overall  stability. 

For  values  of  C  >  0,  corresponding  to  a  negative  sweep  angle,  the  cross¬ 
derivative  term  is  discretized  as: 

c  *1,1+1  ~  *1,1  '  *l-i,l+i  *  *l-i ,J 

j  Ax  Ay 

For  values  of  C  <  0,  corresponding  to  a  positive  sweep  angle,  the  following 
discretization  is  used: 


Ax  Ay 


The  schemes  that  are  presented  for  uniform  mesh  spacing  extend  readily 
to  the  mesh  obtained  from  the  coordinate  transformation.  The  coefficient  of 
the  cross-derivative  is  now 
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For  discretization  of  the  term  c(3£/3x) (3£/3y) (32d>/3C2)  a  centered  scheme  is 
used  at  all  points: 


c . 
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3x 


3? 


i.J  yi+i 


3y 


2*i  +  Vi 
A52 


SOLUTION  ALGORITHM 


The  time-accurate  integration  is  obtained  by  using  an  Alternate  Direction 
Implicit  (ADI)  scheme  that  breaks  the  three-dimensional  problem  into  three  one¬ 
dimensional  problems  in  each  coordinate  direction.  The  advantage  of  this 
scheme  is  its  inherent  stability,  at  least  in  the  case  of  a  linear  equation, 
regardless  of  the  local  Courant  number.  Indeed,  when  solving  a  complicated 
problem,  in  a  mesh  where  the  cell  sizes  may  vary  by  one  or  more  orders  of  mag¬ 
nitudes,  it  would  be  very  time-consuming  to  limit  the  time  step  to  satisfy  the 
Courant  number  associated  with  the  smallest  cell. 

However,  since  the  equation  being  solved  is  nonlinear,  there  is  a  prac¬ 
tical  limitation  which  can  be  associated  with  vortex  shedding  in  lifting 
cases  or  configurations  with  shock  motion.  This  means  that  the  cell  sizes 
must  never  be  so  small  on  the  airfoil  surface  and  near  the  trailing  edge 
that  the  allowable  time  step  for  maintaining  stability  is  unnecessarily 
limited . 


A  Crank-Nicholson  averaging  between  the  levels  n  and  n  +  1  is  used 
since  it  can  be  shown  on  the  linearized  equation  that  a  stable  scheme  results. 

The  three  steps  of  the  ADI-Crank-Nicholson  algorithm  are  as  follows: 


Step  1 
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It  should  be  noted  that  the  underlined  terms  are  treated  explicitly.  However, 
an  implicit  scheme  can  easily  be  devised  based  on  switching  from  a  centered 
difference  approximation  to  32<J>/3£2  when  (c3£/3x)  (3£/3y)  >  0,  and  to  an 
upwind  difference  approximation  when  (c3£/3x) (3£/3y)  <  0.  Test  results  for  a 
swept  tip  showed  very  little  difference  between  the  implicit  treatment  and  the 
explicit  scheme  given  previously. 


Step  2 
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Step  3 
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After  these  equations  are  discretized  according  to  the  method  discussed  in  the 
Finite  Difference  Scheme  section,  the  algebraic  system  is  inverted  by  using  a 
tridiagonal  or  quadradiagonal  direct  solver.  Particular  attention  is  given, 
when  defining  the  finite  difference  analogues,  to  ensure  that  the  main  diagonal 
term  could  be  chosen  as  pivot  in  the  elimination  process.  All  the  terms,  which 
are  treated  implicitly,  contribute  to  the  main  diagonal  with  the  same  sign. 


Each  complete  time  step  requires  approximately  2.5  sec  of  CPU  time  of  the 
CDC  7600  computer.  A  rectangular  blade  computation  requires  approximately 
half  an  hour  of  total  run  time.  For  swept  tips,  where  there  is  a  more  severe 
time-step  limitation,  the  total  run  time  is  an  hour.  The  corresponding  mesh 
is  composed  of  approximately  35,000  nodes. 


RESULTS 


Some  three-dimensional  steady  (hover)  flows  are  simulated  for  three  blade 
geometries:  a  rectangular  blade,  a  swept-tip  blade,  and  a  blade  combining  a 
swept  and  parabolic  tip  (fig.  1).  The  pressure  distributions  are  presented 
for  three  sections  of  the  blades  in  figures  2(a-c)  for  blade  A,  figures  2(d-f) 
for  blade  B,  and  figures  2(g-i)  for  blade  C.  As  can  be  seen,  the  effect  of 
sweep  is  favorable  inboard.  The  shock  waves  either  are  weakened  or  disappear 
on  blade  C.  Close  to  the  tip,  however,  the  opposite  trend  seems  to  occur,  with 
blade  C  experiencing  the  largest  supersonic  pocket.  The  global  effect  is 
favorable  for  the  swept  tips  in  hover. 

Three-dimensional  unsteady  flows  past  a  rectangular  blade  of  aspect  ratio 
/R  *  7,  have  been  computed  at  Mach  numbers  M  =  0.6  and  advance  ratios  (y)  of 
0.45,  0.5,  and  0.55.  The  blade  has  no  twist  and  is  equipped  with  a  symmetric 
NACA  00XX  profile  of  varying  thickness  along  the  span.  ONERA  experimental 
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data  for  the  same  rotors  and  test  conditions  are  available  for  comparison 
(ref.  4).  Figure  3  shows  the  radial  stations  for  the  experimental  pressure 
measurements.  The  corresponding  results  are  shown  in  figures  4(a-c)  for  the 
azimuth  of  60°  and  in  figures  4(d-f)  for  the  azimuth  of  120°  at  the  lowest 
advance  ratio.  For  the  advance  ratio  of  y  =  0.5,  the  results  are  presented 
in  figures  4(g-l).  Figures  4(m-r)  show  results  for  y  =  0.55  for  the  same 
two  azimuth  angles  of  60°  and  120°.  Also  plotted  in  these  figures  are  the 
quasi-steady  results,  which  correspond  to  32<j>/3t3x  =  0.  As  can  be  seen,  the 
unsteady  results  agree  better  with  the  experimental  results,  indicating  a  non- 
negligible  unsteady  term  324>/3t3x.  Furthermore,  a  comparison  of  the  quasi¬ 
steady  solutions  at  azimuth  angles  of  60°  and  120°  exhibits  the  influence  of 
the  cross-derivative  32tj>/3x3y,  which  increases  toward  the  tip. 


THREED  CODE 


THREED  has  been  coded  in  FORTRAN  by  using  only  standard  statements.  In 
its  present  form  it  is  adapted  to  the  CDC  7600  computer  of  the  Ames  Research 
Center,  NASA,  Moffett  Field,  California.  The  Small  Core  Memory  (SCM)  length 
is  27,257  decimal  words  and  the  Large  Core  Memory  (LCM)  length  is  131,072  deci¬ 
mal  words.  THREED  is  divided  into  one  main  program  and  four  subroutines: 

SUBROUTINE  MESH  defines  the  mesh  and  computes  the  metric  coefficients 
SUBROUTINE  SLOPE  computes  the  slope  of  the  blade  at  each  point 
SUBROUTINE  POT  integrates  the  potential  equation 
SUBROUTINE  CP  computes  the  pressure  coefficient  on  the  blade 

The  data  as  they  are  read  in  and  printed  out  are  shown  in  appendix  A.  The 
values  shown  correspond  to  the  results  plotted  in  figure  4(g-l).  A  short 
explanation  of  the  parameters  as  well  as  the  notation  in  THREED  follows: 

i)  HM  =  0. 6  Mach  number 

ALPHO  =  0  mean  incidence  of  the  sinusoidal  motion,  deg 

DALPH  =  0  amplitude  of  the  incidence  variation,  deg 

IROT  =  1  rotating  blade  case 

IROT  =  0  fixed  blade  or  wing  case 

AV  =  0.5  advance  ratio 

GM  =  1.4  ratio  of  specific  heats 

ii)  NSTP  =  601  number  of  time  steps 

ITER  =  400  number  of  iteration  steps 

DTN  =  0.0001  minimum  time-step  size  in  the  relaxation  process  (rad) 

DTX  =  0.01  maximum  time-step  size  in  the  relaxation  process  (rad) 

NMOD  =  8  number  of  elements  in  the  time-step  sequence  based  on 

DTN  and  DTX 

TI  =  -1.5708  initial  time  (rad) 

NPR  =  100  time  step  at  which  results  are  printed 

iii)  YC  “  0.9  location  of  the  kink  or  a  special  span  location 

YD  *  1  tip  of  the  blade 


KSG  =  1 


indicates  that  the  blade  is  symmetric  with  respect  to 
z  =  0 

KSG  =>  0  indicates  that  the  blade  is  not  symmetric 

DEL  =  0.12  thickness  of  the  basic  profile  as  defined  subsequently 

AR  *  7  aspect  ratio  of  the  blade 

iv)  JYSM  *  0  indicates  that  a  strip-theory  condition  is  used  at  the 

root 

JYSM  =  1  indicates  that  a  symmetry  condition  is  used  at  the  root 

KSYM  =  1  there  is  a  lower-upper  symmetry 

KSYM  =  0  there  is  no  lower-upper  symmetry 

KGRAD  =  1  a  Neuman  boundary  condition  is  used 

KGRAD  =0  a  Dirichlet  boundary  condition  is  used 

v)  IM  =  64  number  of  mesh  points  in  the  £  direction  (<64) 

JM  =  32  number  of  mesh  points  in  the  n  direction  (<32) 

KM  =  32  number  of  mesh  points  in  the  C  direction  (<32) 

vi)  IA  =  18  leading-edge  index  (IA  <  IF) 

IF  =  48  trailing-edge  index  (IF  <  64) 

JC  *  11  kink  station  index  (JC  <  JD) 

JD  =  21  tip  station  (JD  <  32) 

KO  ■  16  lower-surface  index 

KU  =  17  upper-surface  index  (KU  -  KO  +  1  <  32) 

vii)  XN  =  -8  location  of  most  upstream  surface 

XX  «  6  location  of  most  downstream  surface 

YN  =  0.5  location  of  innermost  surface 

YX  ■  1.5849  location  of  outermost  surface 

ZN  =  -3  location  of  most  bottom  surface 

ZX  =  3  location  of  most  tip  surface 

viii)  The  basic  airfoil  is  defined  by  two  sets  of  values  representing  the 
abscissas  and  the  ordinates  of  points  on  the  profile.  The  maximum  number  of 
points  is  101.  The  points  are  distributed  in  sequential  order  around  the  air¬ 
foil,  starting  from  the  trailing  edge  describing  the  upper  surface,  then  the 
lower  surface,  and  then  back  to  the  trailing  edge. 

NP  =  101  total  number  of  points  (<101) 

(NP  +  l)/2  must  correspond  to  the  leading  edge 

When  the  geometry  changes  rapidly  it  is  preferable  to  concentrate  the 
points  near  the  leading  and  trailing  edges  in  order  to  ensure  the  best  pos¬ 
sible  accuracy  for  linear  interpolation.  The  coordinate  profile  at  the  mesh 
location  N  is 


XP  (N) 
ZP  (N) 


abscissa 

ordinate 


of  point  N 


The  planform  geometry  of  the  blade  is  defined  by  piecewise-analytic  formulas 
in  the  subroutine  MESH.  An  example  is  given  in  appendix  B  for  a  swept  tip. 
The  functions  XA(J)  and  XF(J)  are  defined  in  the  loop  starting  with 
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77 


DO  6  J  =  1,  JM 


and  ending  with 

87  6  CONTINUE 

as  shown  in  the  box. 


CONCLUSIONS 


A  finite  difference  code  for  predicting  the  high-speed  flow  over  an 
advancing  helicopter  rotor  is  presented.  The  code  solves  the  low-frequency 
transonic  small  disturbance  equation  and  is  suitable  for  modeling  the  effects 
of  three-dimensional  advancing  blade  unsteadiness.  This  work  was  inspired  by 
a  similar  method  developed  by  F.  X.  Caradonna  (ref.  5).  However,  the  computer 
code  THREED  incorporates  some  important  new  features,  especially  tne  capabil¬ 
ity  for  treating  nonrectangular  blade  tips.  Computed  results  show  good 
agreement  with  experimental  blade  pressure  data  and  illustrate  some  of  the 
effects  of  varying  the  rotor  planform.  The  flow  unsteadiness  is  shown  to  be 
an  indispensible  part  of  a  transonic  solution.  It  is  also  shown  that  close  to 
the  tip  at  high  advance  ratio,  cross-flow  effects  can  significantly  affect  the 
solution. 


Ames  Research  Center 

National  Aeronautics  and  Space  Administration 
and 

Aeromechanics  Laboratory 

AVRADCOM  Research  and  Technology  Laboratories 

Moffett  Field,  Calif.  94035,  April  10,  1980 
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APPENDIX  A 


SAMPLE  OF  DATA  AS  READ  IN  AND  PRINTED  OUT  OF  CDC  7600  COMPUTER 


MACH 

NUMBER 

.6 


MEAN  INCIDENCE  ROTATION  ADVANCE  HFAT 
INCIDENCE  VARIATION  Y/N  RATIO  RATIO 

0.  0.  1  .5  1.4 


NO. TIME  NO. STEPS  MIN  STEP  MAX  STEP 
STEPS  RELAXATION  RELAX  RELAX 

601  40  .0001  .01 


RELAXATION  INITIAL 
CYCLE  TIME 

8  -1.5708 


FINAL 

TIME 

1.5708 


SPECIAL  TIP 

SPAN  LOCATION 
LOCATION 

.9  1. 


GEOM.  BASIC  ASPECT 

SYMMETRY  THICKNESS  RATIO 
Y/N 

1  .12  7. 


ROOT  SYM  UP-LO 
CONDITION  SYMMETRY 
Y/N  Y/N 

0  1 


LATERAL 

GRADIENT 

Y/N 

1 


NO  •  X 

NO.Y 

NO.Z 

MESH 

MESH 

MESH 

64 

32 

32 

leading 

EDGE 

18 

TRAILING 

edge 

48 

SPECIAL 
SPAN  NO. 

11 

TIP 

NO. 

21 

LOWER 

SURFACE  NO. 
16 

UPPER 

S .  NO. 
17 

MIN.X 

SURFACE 

-8.0 

MAX  .  X 
SURFACE 
6.0 

MIN.  Y 
SURFACE 

.5 

MAX  .  Y 
SURFACE 
1.5849 

MIN.Z 

SURFACE 

-3. 

MAX  .  Z 
SURFACE 

3. 

IMPR. 

STEPS 

100 
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ABSCISSA  naniMATh 


BASH 
NO.  PI 
101 


:  ppopii.f 

)  I  N  TS  IMDFX 


1 

1.0000 

2 

.9990 

3 

.9961 

4 

.9911 

5 

.9843 

6 

.9755 

7 

.9649 

8 

.9624 

9 

.9382 

10 

.9222 

11 

.9045 

12 

.8853 

13 

.  R64  6 

14 

.8423 

15 

.  R 1  R7 

16 

.7939 

17 

.7679 

18 

.7409 

19 

.7129 

20 

.6841 

21 

.6545 

22 

.6243 

23 

.5937 

24 

.5627 

25 

.5314 

26 

.5000 

27 

.4686 

28 

.4373 

29 

.4063 

30 

.3757 

31 

.3455 

32 

.3159 

33 

.2871 

34 

.2591 

35 

.2321 

36 

.2061 

37 

.1813 

38 

.1577 

39 

.1355 

40 

.1147 

41 

.0955 

42 

.0778 

.0013 

.001^ 

•  00 1  R 
.OOPS 
.0034 
.0046 
.0061 
.0077 
.0096 
.0117 
.0139 
.0163 

•  0 1  RR 
.0214 
.0241 
.  0269 
.0297 
.0326 
.03  54 

.03R2 

.0409 

.0436 

.0461 

.0486 

.0509 

.0529 

.0548 

.0564 

.0578 

.0588 

.0596 

.0600 

.0600 

.0596 

.0589 

.0577 

.0562 

.0542 

.0519 

.0491 

.0460 

.0426 


BASIC  PROFILE 
NO.  POINTS 
101 


INDEX  ABSCISSA  ORDINATE 


4  3 

.0418 

.0389 

44 

.0474 

.0348 

4  6 

.0381 

.0308 

44 

.0248 

.0289 

47 

.0187 

.0211 

48 

.0089 

.0161 

49 

.0039 

.0109 

80 

.0010 

.0088 

81 

.0000 

.0000 

8? 

.0010 

-.0088 

53 

.0039 

-.0109 

84 

.0089 

-.0161 

88 

.0187 

-.0211 

84 

.0248 

-.0289 

87 

.0381 

-.0308 

88 

.0474 

-.0348 

89 

.0418 

-.0389 

40 

.0778 

-.0426 

41 

.0988 

-  .  0460 

42 

.1147 

-.0491 

43 

.1388 

-.0819 

44 

.1877 

-.0842 

48 

.1813 

-.0862 

44 

.2041 

-.0877 

47 

.2321 

-.0880 

48 

.2891 

-.0896 

49 

.2871 

- .0600 

70 

.3189 

-.0600 

71 

.3485 

-.0896 

72 

.3787 

-.0888 

73 

.4063 

-.0878 

74 

.4373 

-.0564 

78 

.4686 

-.0548 

74 

.  8000 

-.0829 

77 

.  8314 

-.0509 

78 

.8627 

-.0486 

79 

.8937 

-.0461 

80 

.6243 

-.0436 

81 

.6845 

-.040° 

82 

.6841 

-.0382 

8  3 

.7129 

-.0384 

84 

.7409 

-.0328 

88 

.  7679 

-.0297 

84 

.7939 

-.0269 

87 

.8  187 

-.0241 

88 

.8423 

-.0214 

16 


BASIC  PROFILE 
NO.  POINTS 
101 


! 


4 


INDEX 

ABSCISSA 

ORDINATE 

89 

.8645 

-.0188 

90 

.8853 

-.0163 

91 

.9045 

-.0139 

92 

.9222 

-.0117 

93 

.9382 

-.0096 

94 

•  9524 

-.0077 

95 

.9649 

-.0061 

96 

.9755 

-.0046 

97 

.9843 

-.0034 

98 

.9911 

-.0025 

99 

.9961 

-.0018 

100 

.9990 

-.0014 

101 

1.0000 

-.0013 
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APPENDIX  B 


SUBROUTINE  MESH 


1)1 1  *=  (  Y  I  J-Y  IC  )  /  (  v  1  O- Y  IC  ) 

Y  (  J  )  =YC+ni)M*  (  YD-YC  ) 
go  Tn  a 

2  CONTINUE 

OI)M=-.S*DI  *  {  Y  I  J-YIC  )  /  Y  IC 
oiiM=cns  (  'him  ) 

Y  (  J  )  =  YN'+nilM*  (  YC-YN  ) 

GO  TO  A 

3  CONTINUE 

.  S#P  !#  I  Y  I  J-Y  1 1)  >/(  I  .  -Y  I  0  ) 

i)iim=cos  ( nu  m  > 

Y  (  J  )  =YX+OUM*  ( YO-YX ) 

A  CONTINUE 

I F ( JSYM.NF.  I  )  Y  (  1  )  =  Y  (  2  ) 

HP  T  T F  ( (S  ,  1000  ) 

00  6  J=1,JM 

WR I TF ( A, 1001  )  J,YI(J),Y(J) 
ft  CONTINUE 

C^sssplanfORM  EQUATION  XA(J),XF(J) 
00  6  J=  1  *  JN 

Y  J  =  Y  (  J  ) 

X A ( J ) =0 . 

X  F  (  J  )  =  I  . 

IF(J.LF.JC)  GU  TO  6 
XA(  J)=3»ft317ft7*(  YJ-YC  ) 

XF( J)=1.+XA( J) 

IF(J.LE.JO)  GU  TO  6 
XA ( J )  =  XA ( J-i  ) 

X F ( J ) =  XF ( J- 1 » 

A  CONTINUF _ 

X  A  (  1  )  =  X  A  I  2  ) 

XF ( 1 )  =  XF ( 2 ) 

WR I TF (ft, 1002  > 

no  7  J=1,JM 

WR I T  F ( ft , 1 00  I )  J , X A ( J  )  ,XF(J) 

7  CONTINUE 
xi t=-i.-oxi 
on  «  i  =  l ,  i  m 
xi i  =  x i i  +  nxi 
xi  <  n  =  xi  i 

ft  CONTINUE 
1=0 
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BLADE  A 


BLADE  A 


(a)  Station  J  =  11. 


(b)  Station  J  =  16. 


(d)  Swept  tip,  station  J  = 


Figure  2.-  Pressure  distribution  computations;  tip  in  hover 


n  =  0.45 . 0.50,  and  0.55 


RECTANGULAR  BLADE 


EXPERIMENTAL  DATA  SECTION  NO.  (|)(3) 


Figure  3.-  Three-dimensional  unsteady  problem  for  forward  flight 


O  EXPERIMENT  O  EXPERIMENT 

—  THEORY  -  THEORY 


(a)  Station  11.  (b)  Station  16. 

O  EXPERIMENT 


(c)  Station  21.  (d)  Station  11. 


Figure  4.-  Pressure  distribution  computations.  Symbols  denote  data  from 

reference  4. 


O  EXPERIMENT 
—  THEORY 


C  =  120 
y  =  0.95 
H  -  0.45 


v  -  120 
V  =  1 
H  »  0.45 


o  EXPERIMENT 


x 


(j)  Station  11. 


O  EXPERIMENT 
—  THEORY 


(k)  Station  16. 


Figure  4.-  Continued. 
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